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ABSTRACT 

This  publication  reports  on  a  compact  presentation  of  a  one-dimensional  model  for  the 
description  of  two-phase  flows  with  heat  exchange.  This  model  was  proposed  earlier  in  the 
form  of  a  system  of  equations,  containing  the  heat  conduction  and  entropy  conservation 
equations  for  each  phase.  In  this  report  the  system  is  formulated  in  a  form  that  is  consistent 
with  the  formulation  of  the  preceding  report  of  the  author.  The  present  system  of  equations 
uniformly  describes  the  evolution  of  averaged  variables  and  interphase  fluctuation 
parameters.  The  correctness  of  the  model  has  been  established  for  convex  equations  of  state 
and  for  proper  compaction  dynamics  kinetics. 
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A  Thermodynamically  Complete  Model  for  One- 
Dimensional  Two-Phase  Flows  With  Heat  Exchange 


Executive  Summary 

Multi-phase  flows  are  common  in  modern  civilian  and  military  environments.  For 
example,  multi-phase  modelling  is  an  inevitable  procedure  required  for  the  design  of 
modern  internal  combustion  engines  and  when  analysing  internal  ballistics  problems. 
In  the  last  decade  the  area  of  multi-phase  modelling  has  rapidly  expanded  due  to  the 
development  of  novel  warheads  dealing  with  heterogeneous  energetic  mixtures,  the 
emergence  of  multi-phase  issues  in  terrorism  threats  such  as  contamination  bombs, 
and  an  extensive  use  of  multi-phase  media  and  porous  materials  for  blast  mitigation. 

Thus,  two-phase  modelling  is  an  urgent  requirement  for  the  analysis  of  novel  weapons 
and  target  responses.  However,  the  nomenclature  of  existing  two-phase  models 
requires  the  use  of  a  specific  software  framework  that  deals  with  a  separate  description 
of  the  two  different  phases  complemented  with  the  pre-selection  of  exchange  terms 
and  connecting  governing  equations. 

The  present  report  deals  with  a  model  presentation  based  on  a  description  of  a  two- 
phase  mixture  as  an  averaged  medium,  using  conservation  laws  for  parameters  of  the 
mixture  complemented  by  governing  equations  for  the  disequilibrium  parameters 
(interphase  fluctuation  parameters).  This  formulation  is  believed  to  be  most  suitable  for 
the  implementation  of  a  model  within  the  existing  hydrocode  frameworks  via  the  user- 
defined  subroutine  apparatus. 

The  present  model  describes  behaviour  of  a  two-phase  heat  conductive  isotropic 
mixture  with  the  ability  to  exchange  mass,  momentum,  and  energy  between  phases. 
The  heat  conduction  is  described  within  the  extended  thermodynamics  approach 
(hyperbolic  heat  conduction)  that  makes  the  final  system  of  equations  of  hyperbolic 
type.  Fundamentals  of  the  model  have  been  published  in  [1],  which  are  slightly 
modified  in  the  present  report  and  reduced  to  a  compact  form  consistent  with  the 
formulation  [2],  Conditions  of  hyperbolicity  of  the  system  of  equations  comprising  the 
model  have  been  established  in  the  report,  which  are  linked  with  the  natural 
conditions  of  convexity  for  the  internal  energy  dependencies  for  each  phase.  The 
correctness  of  the  model  has  also  been  established  by  restrictions  on  the  compaction 
dynamics  kinetics. 
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1.  Introduction 


Multi-phase  media  are  widely  involved  in  modern  civilian  and  military  technology 
environments.  Examples  of  the  media  are  reactive  mixtures  in  high-speed  engines,  a 
variety  of  energetic  and  filling  constituents  employed  in  modern  weapons,  gas-liquid 
aerosols  and  porous  buffer  materials  used  for  blast  mitigation.  In  particular,  the 
heterogeneous  energetic  mixtures  are  widely  used  for  arranging  a  tailored  energy  release 
in  munitions  such  as  delayed  reaction/ initiation,  afterburn,  etc.  Simulation  of  the 
behaviour  of  such  materials  requires  a  description  with  models  which  are  capable  of 
calculating  multi-phase  flows.  The  two-phase  models,  consisting  of  the  conservation  laws 
for  mass,  momentum  and  energy  for  each  phase  are  very  popular.  In  that  case,  the 
conservation  laws  are  formulated  for  partial  characteristics  (the  mass-related 
characteristics  for  each  phase  with  respect  to  the  mixture  volume,  containing  several 
phases)  and  interconnected  by  the  exchange  terms.  The  system  of  equations  for  both 
phases  is  complemented  by  a  tabulated  governing  equation  for  the  volume  concentration 
of  a  phase  (the  compaction  dynamics  equation).  This  approach  has  certain  drawbacks:  i)  it 
involves  non-conservative  terms  associated  with  so-called  'interfaciaT  (or  configurational) 
pressure  (e.g.,  see  the  frequently  cited  Baer-Nunziato  (BN)  publication  [3]  for  a  model  that 
was  first  formulated  in  [4]);  and  ii)  it  is  impossible  to  implement  the  model  in  a  hydrocode 
with  the  user-defined  subroutine  interface  because  the  description  involves  descriptions 
for  each  phase  and  it  is  in  contradiction  with  the  typical  hydrocode  structure.  Many 
attempts  to  consider  a  multi-phase  medium  as  an  averaged  one  have  been  made,  including 
a  classic  monograph  by  Truesdell  [5].  However,  a  closed  thermodynamic  formulation, 
resulting  in  an  efficient  practical  realisation,  had  not  been  proposed  at  that  time.  A  variety 
of  models  have  been  recently  developed  in  several  papers  [3,  4,  6,  7,  8],  However,  they  are 
not  formulated  as  conservation  laws  and  these  non-conservative  formulations  complicate 
analysis  of  thermodynamical  correctness  of  the  models. 

The  present  work  employs  an  approach,  enabling  us  to  derive  equations  in  the  form  of 
conservation  laws.  One  of  the  first  realizations  of  this  approach  has  been  published  in  [9]. 
This  allowed  one  to  formulate  equations  for  the  averaged  parameters  and  mass 
concentration,  so  the  effective  averaging  parameter,  which  is  involved  in  the 
consideration,  is  the  mass  concentration  of  a  phase.  Realization  of  this  approach  as  a 
computer  code  [9]  resulted  in  significant  numerical  difficulties  associated  with  the  mass 
exchange  in  the  areas  of  high  mass  concentration  due  to  actual  involvement  of  partial 
characteristics  in  the  calculation  process.  A  formulation  [7]  should  be  mentioned  which 
modified  the  BN  model  into  a  form  that  combines  three  conservation  laws  for  the  mixture 
with  four  independent  equations  out  of  the  original  BN  model.  Replacement  of  the  partial 
characteristics  by  specific  characteristics  of  the  phases  ('real'  densities,  pressures,  etc) 
resulted  in  the  parametrization  (see,  e.g.,  [10])  that  claims  that  the  mass  concentration  and 
density  are  not  the  only  parameters  involved  in  the  complete  description  of  the  phases' 
state.  This  parametrization  dictates  that  characterisation  of  the  phases'  state  has  to  involve 
both  the  mass  and  volume  fractions.  The  models  employing  this  parametrization  resulted 
in  a  successful  model  such  as  [11].  Reduction  of  the  model  for  the  one-dimensional  multi¬ 
phase  flows  to  a  compact  form  has  been  conducted  in  [2], 
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The  model  [2,  11]  considers  a  two-phase  flow  without  heat  exchange,  by  neglecting  the 
temperature  and  entropy  differences  between  the  phases.  However,  the  heat  exchange 
appears  to  be  critical  in  mixtures  that  contain  phases  which  are  either  highly  energetic  or 
have  highly  contrasting  properties  and  very  different  thermo-physical  characteristics.  A 
model  that  takes  into  account  the  heat  exchange  on  the  basis  of  hyperbolic  heat 
conductivity  was  formulated  in  [1],  The  equation  derivation  in  the  present  report  modifies 
the  model  in  order  to  preserve  the  characteristic  structure  of  the  original  systems  when 
reducing  it  to  the  two-phase  case.  The  outcome  suggests  a  reduced  form  of  the  model  for 
the  one-dimensional  case,  and  proposes  a  set  of  parameters  allowing  one  to  derive  a 
compact  form  of  the  model  in  this  case.  This  formulation  helps  us:  i)  to  provide  a 
thermodynamically  correct  model;  ii)  to  solve  a  Riemann  problem  for  a  two-phase  mixture 
with  the  help  of  a  single-phase  Riemann  solver  applied  to  each  phase;  and  iii)  to  formulate 
jump  conditions  in  a  compact  form.  This  formulation  could  also  be  convenient  for 
generalisation  of  the  model  to  multi-phase  case. 


2.  Basic  single-phase  model 


To  start,  we  recall  a  heat  conduction  single-phase  model  [12]  for  the  case  of  isotropic 
material  with  a  spherical  stress  tensor  presented  by  pressure  only.  The  model  is  based  on 
the  principles  of  extended  thermodynamics  [13]  and  uses  a  generalisation  of  Fourier's  heat 
conduction  law  to  a  non-stationary  heat  conduction  hyperbolic  equation.  We  slightly 
generalise  the  model  by  including  the  right-hand  side  source  terms  that  characterise  the 
external  mass,  momentum,  and  heat  production  (this  generalisation  will  be  useful  further 
on,  when  the  system  is  used  for  the  development  of  a  two-phase  model).  The  combined 
system  of  equations  for  the  description  of  mechanical  and  thermal  response  of  a  material 
takes  the  following  form  [12]: 


dp  dpu 

—  +  -A—  =  m 

8t  dx 


0  ’ 


dpu  d{pu 2  +  p ) 


dt 


dx 


dpq  d{pqu  +T) 
dt  dx 

dp(e  +  u2  /l)  d\pu[e  +  u2 /l) 


l0  5 

pq 

T 


dt 


dx 


+  pu  +  TJ  7 
—  1 0  , 


(1) 


here  mo  is  the  mass  production  rate,  m  is  the  momentum  production  rate  and  lo  is  the 
energy  production  rate.  Dependent  thermodynamic  variables  are  calculated  as 

p  =  p2eP ,  T-eSr  J  =  eq.  (2) 

They  are  associated  with  the  following  thermodynamic  identity 
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T  ds  -  de  +  p  dV  -  j  dq 

where  internal  energy  e  =  e(p,  s,  q)  is  a  given  function  of  the  material  state.  Independent 
variables  of  the  system  (1)  are  p,  s,  q,  and  u,  which  are  specific  density,  specific  entropy,  a 
thermal  impulse  variable  and  particle  velocity,  respectively.  The  function  r  determines  the 
characteristic  time  of  'thermal  relaxation'  (see  [12]).  The  dependent  variables  are  pressure 
p,  temperature  T,  and  an  entropy  flux  variable  /,  which  are  calculated  from  the  internal 
energy  as  stated  in  (2).  It  should  be  noted  that  the  origin  of  the  variable  q  comes  from  the 
heat  flux  entering  Fourier's  law,  but  in  the  present  form  it  has  quite  a  vague  physical 
meaning  with  units  in  (s-m2-grad)/kg.  In  fact,  its  physical  meaning  is  that  /  =  eq  (/  is 
measured  in  W/ (m2-grad)),  where  the  actual  heat  flux  is  presented  by  the  term  J-T. 

Expanding  the  last  equation  of  the  system  (1)  with  the  use  of  the  mass  and  momentum 
conservation  laws,  we  have 


de  du  dTJ 

p —  +  p  —  + - 

dt  8x  8x 


=  h 


here  d/dt  =  d/dt  +  u-d/dx  is  the  particle  derivative.  Using  the  calculation  rules  for 
dependent  variables  and  the  heat  flux  equation  of  the  system  (1),  we  can  derive  the 
entropy  production  equation: 


dps  dpsu  +  J 
dt  dx 


-m0/j'  -un0 


+  <pd), 


(3) 


where  cp  =  pq/r  and  p '  =  p  -  u 2/  2  when 

p-e+p/p-Ts-qj  (4) 

is  the  chemical  potential.  The  equation  (3)  shows  that  J  is  really  the  entropy  flux  and  the 
entropy  production  is  non-negative  (in  the  absence  of  external  mass  and  momentum 
production)  when  q-J>  0. 

For  the  derivation  of  the  system  of  equations  in  the  form  of  conservation  laws,  one  more 
consequence  of  the  system  (1)  will  be  useful  when  we  deal  with  the  two-phase  case.  We 
note  that  a  consequence  of  the  momentum  equation  could  be  rewritten  as 

du  du2  2  1  dp  nn-umn 

—  + - —  + - —  =  — - -  .  (5) 

dt  dx  p  dx  p 

This  equation  has  been  used  in  [1,  9,  11]  in  order  to  derive  an  equation  for  the  velocity 
interphase  fluctuations  in  a  two-phase  single-temperature  non-heat-conducting  medium. 
This  was  possible  because  the  non-conservative  term  in  the  left-hand  side  of  (5)  could  be 
reduced  to  a  gradient  of  the  chemical  potential  because  the  temperature  variation  between 
the  phases  was  neglected.  In  the  present  case,  using  the  derivation  rules  (2)  for  the 
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dependent  variables,  the  differential  of  the  chemical  potential  (4)  can  be  assessed  as 
follows 


dp '  =  de  + 


dp 

P 


dp-Tds  -  sdT  -  qdJ  -  Jdq 
P~ 


—  -  sdT  -  qdJ  . 
P 


(6) 


From  (6),  equation  (5)  can  be  rewritten  in  the  following  form 


du  du2  jl  dp  dT  8J  n0-um0 

—  + - — +  — +  5 - +  q —  =  — - - 

dt  dx  dx  dx  dx  p 


(7) 


It  is  clearly  seen  from  (7)  that  for  a  single-temperature  medium  (when  temperature  and 
entropy  are  not  changing  between  the  phases)  in  the  absence  of  heat  conduction,  the 
velocity  difference  allowed  one  to  keep  the  chemical  potential  in  the  divergent  form  and 
eliminate  the  temperature  and  entropy  flux  gradients  from  the  velocity  difference 
equation.  For  the  heat  conductive  medium  with  varying  temperature  between  the  phases 
the  situation  is  more  complex.  However,  it  is  possible  to  derive  another  conservation  law 
using  a  consequence  from  the  heat  conduction  equation  of  the  system  (1)  rewritten  for  q 
and  the  entropy  conservation  law  (3).  Multiplying  the  heat  conduction  equation  by  ps  and 
the  equation  (3)  by  q,  after  their  summation  we  have: 


8 psq  dpsqu 
8t  dx 


8T  8J 

+  s - 1-  q  — 

dx  dx 


q  (j  ,  \  qJ-sT 

~Vo -m0p  -un0  -sTm0)+  (p — - — 


and,  subtracting  this  equation  from  (7),  we  obtain  a  required  conservation  law 
d{u  -  psq )  s(w2/2  -  psqu  +  p) 


dt 

n0  -um0 
P 


dx 

q  (, 

-(do  ~mo P  ~uno 


\  qJ-sT 
-sTm0)-<p — - — . 


(8) 


For  the  sake  of  convenience  we  call  the  difference  u  -  psq  a  chemical  flux  due  to  the 
association  between  the  variable  change  and  the  gradient  of  p'\ 


dpq  d(prju  +  p1) 

dt  dx 


n0  -um0 


q  p  ,  rr  \  qJ-sT 

-  — l/o  ~m0p  -un0-sTm0)-(p - - — , 

p  T  T 


where  prj  =  u  -  psq. 

The  system  of  equations  (1)  complemented  by  the  auxiliary  equations  (3)  and  (8)  and 
closed  by  relationships  (2)  are  all  necessary  to  derive  equations  for  a  two-phase  model. 
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3.  Selection  of  parameters  for  a  two-phase  model 


Let  us  denote  the  averaged  density  of  a  two-phase  medium  by  p  =  m/V,  here  m  is  mass  of 
a  representative  volume  and  V  is  the  volume  occupied  by  the  mass.  Similarly,  we  can 
define  specific  densities  of  the  phases  pi  =  mi/Vi  and  pi  =  1112/V2,  where  m  =  mi  +  m2  and 
V=  Vi  +  V2 .  Multi-phase  theories  usually  deal  with  so-called  partial  densities,  which  relate 
the  phase  masses  to  the  whole  volume  such  as:  p'\=m\/V  and  p' 2=1)12/  V.  These  partial 
characteristics  are  key  ones  because  the  conservation  laws  for  each  phase  are  actually 
formulated  only  for  the  partial  characteristics.  For  the  case  of  a  media  with  phases,  which 
are  capable  of  exchanging  mass,  momentum,  and  heat,  the  conservation  laws  in  the  one¬ 
dimensional  case  take  the  following  form  for  the  first  phase: 


dp[ 
dt 
8p[u  1 
dt 

dp[  q  1 


1  dp\ux 
+  =  m 


dx 


0  ’ 


=  n 


+ 


d{p[ul  +p[) 
dx 

d{p\  qxux  +  Tx ) 


0  ’ 


A 


dt  dx  Tj 

dp[(e[+u2j2)  |  d[p[ux{e\  +  u2j2) 


+  p,Wj  + 


T  J 

1  1 


dt 


dx 


=  /n 


(9) 


and  for  the  second  phase: 


dp '  dp\u , 

-+  =-m 


dt 
dp'2u2 
dt 

dp'2q2 


dx 


0  ’ 


+ 


+ 


=  ~n 


d[p\u\  +  P2) 
dx 

d(p2  q2u2  +T2) 


0  5 


Pi  $2 


(10) 


dt  dx  t2 

dp2(e'2+u22/  2)  d[p2u2(e'2  +  u22/2)+ p'2u2  +T2J2\ 


dt 


dx 


=  -/„ 


Here  mo,  no,  and  lo  are  the  mass,  momentum,  and  heat  exchange  rates,  respectively;  u, 
(i=l,2)  are  velocities  of  the  phases,  p T  and  p'2  are  partial  pressures  within  the  phases,  cji 
(i=l,2)  are  thermal  impulse  variables  for  the  phases,  and  e\  and  e  2  are  specific  internal 
energies.  The  thermodynamic  identity 


T  ds-  de  +  p  dV  -Jdq-de-p  dp/p  2  -Jdq,  (11) 

where  V  is  specific  volume,  being  applied  to  each  of  the  phases,  enables  us  to  calculate 
partial  pressure,  temperature,  and  the  entropy  flux  similar  to  (2): 
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>;=(a'): 


Pi 


i=(p'i)2 


de[ 

to’ 

:  de2 

to’ 


T  -M 
1 

T2=  — 

ds. 


J 1  = 


5e[ 


J  2  ~ 


Se2 

dq2 


(12) 


if  dependencies  of  specific  energy  on  p  ,  s,  and  q  are  given: 

ei  =  e'i(p'i,  si,  qi)  ,  ei- e  p  i,  si,  qi)  .  (13) 

It  should  be  noted  that  the  definition  of  partial  pressure  is  based  on  the  application  of  the 
thermodynamic  identity  (11)  with  respect  to  partial  density.  Thus,  the  traditional  approach 
to  the  calculation  of  two-phase  flows  is  to  solve  numerically  systems  (1)  and  (2),  pre¬ 
selecting  the  exchange  terms  mo,  no  and  lo,  and  tabulating  the  'equations  of  state'  in  the 
form  (13)  (for  the  sake  of  convenience,  we  call  the  relations  such  as  (13)  equations  of  state), 
using  (12)  for  the  calculation  of  pressure,  temperature,  and  the  entropy  flux. 

We  can  complement  the  systems  (9-10)  with  additional  conservation  laws  for  the  entropy 
flux  similar  to  (3) 


dp[sx  dp[sxux  +JX 
dt  dx 

dp'2s2  |  dp'2s2u2  +  J2 
dt  dx 


~ (/<>  - m0p[  ~uxn0  +  (pxJx)  =  Rx  , 

A 

=  /„  +  moPl  U2  n0  (Pl^l  )  =  ^2  ’ 


(14) 


and  for  the  chemical  flux  as  follows  from  (8): 

d{ux- p[sxqx)  |  d(ux  /2-  p[sxqxux  +//,)_ 
dt  dx 


-qr(/o  ~mof*[  -Uin0  -sxTxm0)-<px 


A 


qxJ\  sxTx 


A 


=  N, 


i  ■> 


d(u2~  p'iSiqi)  |  d(u22/2- p'2s2q2u2+ p2)  _ 
dt  dx 

/In  U2m0  (l  t  rri  \ 

=  -  b  (/(,  —  m0p2  —  u2  n0  —  s2T2m0 )  —  cp2 

Pi  T2 


q2J2  s2  T2 


=  n2. 


(15) 


It  should  be  noted  that  the  chemical  potentials  p\  may  involve  in  their  definitions  (that  are 
identical  to  (4))  equally  partial  (p'\  and  p'\)  or  'real'  (p,  and  pi)  characteristics,  because  p,/ pi  = 
P'i  Ip  i  (see  [2]). 

The  procedure  of  averaging,  having  been  employed  in  [7],  involves  the  introduction  of  the 
averaged  density,  pressure,  and  velocity.  We  introduce  as  in  [10]  the  mass  concentration  of 
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the  first  phase  as  c  =  C\  =  mi / m,  then  for  the  second  phase  ci  =  mi/  m  =  1  -  c.  Then  we  can 
calculate  the  partial  densities 


p'i=mi/V=(mi/m)-(m/V)=pc  ,  p'2-p(l-c)  .  (16) 

Introducing  volume  concentration  for  the  first  phase  as  0i  =  0  =  Vi/V  and  for  the  second 
phase  as  02  =  1  -  0  =  Vi/  V,  we  can  recalculate  the  specific  densities  of  the  phases  as  follows 

pi  =  mi/V i=  (mi/ m)  -(in/  V)  {V/Vi)  =  pc/d  ,  pi  =  p(  1  -  c)/(  1  -  0) .  (17) 

It  is  common  (see  [3,  4,  6-9,  11])  to  tabulate  a  governing  equation  for  the  parameter  0  (the 
compaction  dynamics  equation)  evolving  0  along  the  particle  path  of  a  mixture.  Then  from 
(15)  an  evolutionary  equation  for  c  in  the  form  of  the  conservation  law  follows  from  the 
conservation  of  mass  for  the  first  or  second  phase.  Thus,  selection  of  the  parameters 
related  to  density,  namely,  the  averaged  density  and  mass  and  volume  fractions  can  be 
arranged  in  a  natural  way. 

Let  us  calculate  pressure  in  an  averaged  medium.  Firstly,  we  link  the  local  pressures  and 
densities  with  the  partial  ones.  We  consider  an  alternative  presentation  to  (13)  for  the 
equations  of  state  in  the  following  form 

ei  =  ei(  pi,  si,  qi)  ,  e2  =  e2(  p2,  s2,  q2)  .  (18) 


Then,  from  (12)  and  (16)  it  follows  that 


del 
dp[ 

P\  =(l  ~0)Pi  > 


A=(A')2fV  =  M): 

dpi 


1  de1 


6  dp, 


.2  de, 


=  e-pt-^  =  e-Pl 

dp\ 


(19) 


where  the  following  denotations  for  the  local  pressures  are  used: 


Pi  =  P\ 


2  de, 


dpi 


Pi  =  Pi 


2  de2 


dp2 


(20) 


Obviously,  choice  of  the  equation  of  state  in  the  form  e  or  e  does  not  affect  the  calculation 
rules  (12)  for  temperature  T  and  entropy  flux  /. 

Selection  of  the  parameters  for  the  variables  associated  with  momentum,  thermal  impulse, 
and  entropy  is  not  so  easy.  Because  our  objective  is  to  derive  equations  in  the  form  of 
conservation  laws,  it  seems  obvious  to  select  parameters  from  combinations  of  the  basic 
variables  entering  the  conservation  laws  for  the  phases  (9)  and  (10). 


We  start  with  the  momentum  conservation  law  that  gives  for  the  total  momentum  the 
combination  pT  Mi  +  p'2  u2.  Alternative  combinations  resulting  in  conservation  laws  are 
either 
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p'lUi-p'iui  (21) 

or 

ui- p\siqi-U2  + p'isiqi  .  (22) 

The  last  combination  is  in  agreement  with  a  conservation  law  obtained  from  equations 
(15).  From  the  combination  for  the  total  momentum  conservation  it  follows 

pu  -  p'lUl  +  p'lUl 

and  from  (16) 

u  =  c  ui  +  (1  -  c)  ui .  (23) 

Selection  of  parameters  for  the  thermal  impulse  and  specific  entropy  may  follow  from  the 
characterisation  of  the  total  thermal  and  entropy  fluxes 

piqi  +  p'2^2 

p'isi  +  p'  2  S2 

and  alternative  combinations  follow  from  the  fluxes'  differences: 

p'iqi-p'iq2  (24) 

and 

p’lSl  -p'2S2  .  (25) 

A  natural  restriction  follows  from  the  energy  conservation  laws  in  (9)  and  (10),  from  which 
the  total  energy  is  calculated  as  follows 

p'l  (ei  +  ui2/2)  +  p'2(e2  +  «22/2)  .  (26) 

Because  specific  energy  is  an  extensive  variable,  the  averaged  specific  energy  is 
e  -  c  ei  +  (1  -  c)  P2  .  (27) 

Thus,  it  follows  that  (26)  and  (27)  should  give 


7  2  X 

(  2  ^ 

(  2  X 

^  a  1 

u ,  1 

E  +  — 

=  A 

ex  +  — 

+  Pi 

l  2  7 

l  2  J 

l  2  J 
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where  u  is  determined  from  (23).  Then,  using  (16),  we  have: 


E  H - =  c 

2 


.2  A 


e\  +■ 


+  0-d 


.2  3 


e2  +- 


=  e  +  - 


cu: 


+  (l  -  c)i 


Therefore,  using  (23),  we  must  define  a  generalised  internal  energy  E  as  follows 


CU2  +(l-c)«2  -{cul  +(l-c)w,)2  c(l  -c)(ul  -U2)2 
E  =  e  H - - - - —  =  e  H - 


(28) 


It  should  be  noted  that  the  definition  (28)  is  not  dependent  on  how  a  choice  of  variables 
affects  the  velocity  difference.  The  only  restrictions  imposed  were  those  on  the  total 
momentum  (23)  and  on  the  averaged  specific  internal  energy  (27). 

To  choose  a  parameter  for  the  velocity  fluctuation  we  have  two  options  (21)  and  (22).  We 
start  our  analysis  with  the  option  (21),  denoting  a  velocity  difference  parameter  by  w,  so 
that 


pw  -  p  i  Wi  -  p  2«2 

or,  using  (16)  and  (23),  we  have 


u  +  w 


u-w 


(29) 


1_  2c  ’  2  ~  2(l  -c) 

This  definition  reduces  the  generalised  internal  energy  (28)  to  the  following 
F_r.  [(l-2c)i/  +  w]2 

4c(l-c)  ' 

Subtraction  of  the  momentum  equations  from  (9)  and  (10)  along  with  (19)  and  (29)  gives 


|  u~  U- 


Sf  y ouw  + 

Spw  +  _[ _ 

dt 


(l  -  2c)  +  (2  -  4c  +  4c2 )  uw  +  w2  (l  -  2c)] 


4c(l  -c) 


d[0px  -  (l  -  d)p2  ] 


dx 


dx 


=  2 n0  . 


Thus,  this  governing  equation  clearly  shows  that  evolution  of  the  parameter  w  is 
determined  by  the  gradient  of  a  function  which  can  not  be  associated  with  a  potential  of 
the  generalised  internal  energy  E.  For  this  reason  this  parametric  choice  was  considered  to 
be  unsuitable  for  our  purposes. 

Thus,  the  only  viable  option  for  the  selection  of  parameters  characterising  momentum, 
thermal  impulse  and  entropy  are  combinations  for  the  total  momentum,  total  thermal 
impulse  and  total  entropy  flux  complemented  with  the  combinations  (22),  (24),  and  (25). 
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With  the  use  of  (16)  for  the  thermal  impulse  and  entropy  flux,  these  combinations  result  in 
the  following  equations: 

q  =  cql+(l-c)q2,  A  =  cql  -(l- c)q2  , 

S  =  CJj  +  (l  -  c)s2  ,  X  -  CSl  -  (l  ~C)S2  • 


A  generalisation  for  multi-phase  (n-phase)  case  might  take  the  following  form  for  Ai  (an 
analogue  of  the  parameter  A):  Ai/ n  =  Ci  qi  -  q/ n  (i=l,. .  .,n).  It  is  easy  to  see  that  for  the  two- 
phase  case  this  definition  is  reduced  to  the  definition  (30):  A  =  2c  qi  -  [c  qi  +  (1  -  c)  qi]  or  A  = 
c  qi  -  (1  —  c)  qi.  Thus,  the  phase  thermal  impulse  and  entropy  can  be  calculated,  using  the 
present  choice  (30)  that  results  in 


q  +  A  q-A  s+x  s~x 

q'=—'  q>=Wcj'  St=—’  Si=w^r 


(31) 


Using  the  definitions  (31),  (16)  and  (23),  and  the  combination  (22)  that  was  selected  for  our 
derivation,  we  can  write  down  definitions  of  parameters  involving  velocities: 


u  =  cu{ 


+  (l  -c)u2  , 


pv  =  Wj  -u2 


pUq  +  A)(s  +  x) 
4  c 


(q  ~  AXS  ~  x) 

1  -c 


(32) 


For  the  sake  of  convenience  we  denote 


v  =  v  +  — 
4 


(q  +  A)(s  +  x)  (q-X)(s-z) 


I  -  c 


(33) 


then  from  (32)  and  (33)  the  phase  velocities  are 
«i  =  u  +  p  v '  (1  -  c)  ,  ui  =  u  -  pv'  c  . 

Summarising,  the  generalised  energy  can  be  calculated  from  (28)  as  follows 

E  =  e+4-c)(Pv')2 


(34) 


(35) 


where  v ",  which  is  defined  by  (33),  is  a  function  of  v,  q,  s,  c,  A,  and  j.  It  should  be  noted 
that  with  the  present  choice  of  velocity  parameters  (32)  the  variable  v  does  not  have  a 
dimension  of  velocity,  similarly  to  the  definition  of  velocity  difference  in  the  compact 
formulation  of  the  model  [1]  obtained  in  [2], 
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4.  Equation  of  state  for  the  two-phase  model 


Using  the  equations  of  state  (18)  for  the  phases  and  definitions  (17)  and  (31),  we  can  build 
up  the  equation  of  state  for  an  averaged  medium: 

e(p,  c,  Q,  q,  A,  s,x)  =  c  ei(pc/6,  (q  +  A)/  (2c),  (s  +  j)/  (2c))  + 

(36) 

+  (1  -  c)  e2(p(  1  -  c)/(l  -  0),  (q  -  A)/ (2(1  -  c)),  (s  -j)/(2(l  -  c)))  . 

Summarising  the  chosen  set  of  parameters,  the  kinematic  variable  of  the  model  is  averaged 
velocity  u  and  internal  variables  (independent  thermodynamic  parameters)  are  p,  c,  9,  q,  A, 
s,  j,  and  v.  Using  (35)  along  with  (36),  we  can  calculate  derivatives  of  E  with  respect  to 
each  of  the  independent  thermodynamic  parameters.  We  start  with  the  derivative  of  E 
with  respect  to  the  mass  concentration  c: 


Ec  =  e,  -e2  +c 


dgt  P 

dpi  0 


-(l-c) 


de2  p  dex  q  +  A 


-c 


dp2  1  —  6  dql  2 c 


+  (l-c) 


de2  q  -  A 


dq2  (l-c)2 


—  c 


dex  s  +  x  i,  \de0  s-%  (l-2 c)(pv')2  (  , 

^iTT+(1-£UW!+! - +  -c)p 


ds,  2c 


ds2  (l-c) 


—  Cj  e2  +  px 


de. 


dpi 


Pi 


de~, 


de, 


2 

de0 


,  8v' 

dc 


v 


■  +  5 


-  ^2  -  S1  ~  ■  ~z  ~ 

op2  oqx  oq2  os ,  dss 


dex  de2  (l  -  2c){pv') 


■  + 


-c(l-c) 


p2v' 


(q  +  A)(s  +  z)  (g-A)(s-z) 
(l-c)2  . 


Combining  the  result  with  (17),  (20),  (31),  (12)  and  (4),  we  obtain  the  derivative  in  a 
compact  form,  which  may  be  called  a  generalized  chemical  potential: 


A  =  Ec  =  Mi  -  Pi  + 


(l-2c)(/V)2 


-c(l-c)^ 
V  4 


(q  +  A)(p  +  x)  {q  ~  A)(s  -  x) 
(l-c)2  . 


(37) 


It  should  be  mentioned  that  a  similar  inverse  dependence  of  phase  entropies  on  the  mass 
concentration  chosen  in  [14]  proved  to  be  a  suitable  representation  for  the  potential  Ec  to 
be  associated  with  an  affinity  of  the  Gibbs  potentials. 

The  derivative  of  E  with  respect  to  the  specific  density  p  can  be  associated  with  a 
generalised  pressure  P  and  it  is  calculated  as  follow 
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r.  2n  2  ^  C  /.  \  <9<?2  1  -  C 

P  =  P'Ed  =  p  c - +  (l-c) — - - 

'  l  dp,  o  v  ’ dPl  \-e 

=  Op,  +  (i  -  e)p2  +  c( i  -  c)p3  (v')2  • 


+  c(l-c)/?3(v')2  = 


(38) 


The  derivative  of  E  with  respect  to  0  can  be  calculated  with  the  use  of  (17),  (20)  and  (36): 

(39) 


n  =  Ee  =  -c^+ (i  -  c)P^Pl  -d  =  -  PcLgL . 


p 


A  generalised  characteristic  associated  with  the  entropy  flux  /  could  be  obtained  by 
differentiation  of  E  with  respect  to  q: 


„  de,  1  /  \  de, 

J  =  En  =c^—  +  (l-c)  2 


2c 


dq2  2(1  -c) 


1  h  \  2 

+  c(l-c)/rv  — 

oq 


P-j,P-j1+4-c)PP(s-^-^P 

2  2  v  '  4  ^  c  1  —  c 


(40) 


and,  similarly,  generalised  temperature 


T  =  E=c 


de,  1 
8s,  2  c 


■  + 


(i-c) 


de. 


ds2  2(l -c) 


1  h  \  2  ,8v' 

+  c{l-c)pv  — 

ds 


P-TtP-T1+c(l-c)PP^  +  X  q-P 

2  2  4 


V  c 


1  -  c 


(41) 


Derivatives  with  respect  to  A  and  j  can  be  obtained  in  a  similar  way: 

- — bc(l-c)/r 


^  de,  1  /  \  de, 

B  =  E,  =c^ — (l  —  c) — - 

dA  c  d^2  1  -  c 


■\'8v‘ 


dA 


=  -J, -~J2  +  c(l-c)^ 
2  2  4 


2"'rs  +  x  | 


V  c 


1  —  c 


(42) 


and 


A=  E  =cgfLl-(l-c)^  1  +  c(l-c)p2v 


ds,  c 


ds2  1  -c 


,dv' 


=  -T,  -~T2  +  c(l-c)^~ 
2  2  4 


q  +  A  +  q  -  A^ ) 


V  c 


l-c 


(43) 


Temperatures  and  entropy  fluxes  of  the  individual  phases  can  be  calculated  from  (41),  (43) 
and  (40),  (42),  respectively,  as  follows 
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T,=T  +  A-c{l-c)^-^,  T,=T-A+c(\-c)^-^—2 
2  C  2  1  ~C 


(44) 


and 


2  t 


P2v'  s-x 


J,  =J  +  B-c(l-c)^^-^±^,  J,  =J-B  +  c(l-c) 

2  c  2  1  —  c 

Finally,  from  (35)  we  need  to  calculate  the  derivative  of  E  with  respect  to  v: 

Q  =  Ev  =  c(l  —  c) p~v'  ——  =  c(\-c)p2v' . 

dv 


(45) 


(46) 


This  concludes  the  set  of  formulas  for  the  determination  of  the  dependent 
thermodynamical  parameters  of  the  model. 


5.  Conservation  laws  of  the  two-phase  model 


Conservation  laws  will  be  derived  from  the  conservation  laws  (9)  and  (10)  for  partial 
parameters.  We  start  with  the  mass  conservation  law.  From  (16)  and  (34):  p\  +  p\  =  p  and 
p\  Hi  +  p' i  Ui  =  p  u.  Using  these  relations  and  summing  up  the  continuity  equations  in  (9)- 
(10),  we  can  obtain  the  continuity  equation  for  the  averaged  variables: 


8P  !  dpu  _ Q 
dt  dx 


(47) 


Rewriting  the  continuity  equation  in  (9)  with  the  use  of  (17)  and  (34),  we  can  obtain  the 
following  governing  equation  for  the  mass  concentration: 


dpc  d[ 


puc  +  p  c 


(l -cV 


dt 


dx 


=  m 


o  • 


From  (46)  this  equation  can  be  reduced  to  the  following: 


dpc  dipuc  +  Q.) 

■  +  — — - -  =  m 


dt 


dx 


0  ’ 


(48) 


where  Q  =  Ev. 

The  compaction  dynamics  equation  for  the  volume  concentration  is  chosen  in  the  well- 
know  form  of  a  conservation  law  within  a  material  volume  [3,  4,  6-9, 11, 14]: 
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dp6  dpuO  ,  pU 

-t—  +  — —  =  $>  =  -—,  (49) 

dt  dx  r0 


here  O  is  a  function  responsible  for  the  process  of  phase  compaction  and  to  is  the  relaxation 
time  function  governing  the  compaction  kinetics.  It  should  be  mentioned  that  the  kinetics 
chosen  in  (49)  is  quite  natural  because  the  compaction  dynamics  (rate  of  change  of  the 
volume  fraction)  should  be  associated  with  the  pressure  difference  between  the  phases. 
This  association  can  be  easily  seen  from  the  definition  (39)  for  n. 

The  momentum  conservation  laws  (the  second  equations  of  systems  (9)  and  (10))  give  the 
following  momentum  equation  for  the  two-phase  medium: 

dp!L+iPic+p)=(>  (50) 

dt  dx 


Here  the  following  relationships  have  been  used  that  were  derived  from  equations  (34) 
and  (38): 


p'l  Ml  +  p'l  U2- p  u  , 

p'l  (Ml)2  +  p'l  (lll)2  +  p  =  pil2+p  +  p3  C(1  -  c)(v')2  =  p  U2+  P  . 

As  we  discussed  in  Section  3,  auxiliary  conservation  laws  (15)  for  chemical  flux  will  be 
used  for  the  representation  of  the  velocity  difference  between  the  phases.  We  subtract  the 
second  equation  of  (15)  from  the  first  one  and  for  the  first  term  of  the  left-hand  side  we  use 
the  representation  (22)  directly  in  the  form  (32).  The  second  term  of  the  left-hand  side  of 
the  resultant  equation  can  be  calculated  with  the  help  of  the  following  relationships: 


(Ml)2  -  (M2)2  =  (mi  -  Ui) (ui  +  M2)  =  2 puv'  +  (1  -  2c)  ( p  v')2  , 


t  Pu 

PU  V  =  puv  +  -j- 


(q  +  2)(s  +  x)  (q~  A)(s  -  z) 


1  -c 


and 


p2s2q2u2  -  plsxqlux  = 


P 


k  -  z)U 1  -  h(„  _  pv'C)-k±Akid(u + pv'{\-c)) 


1  -c 


pu 


(s - x)(q -a)  {s  +  x){q  + 


1  -c 


-c(i-c) 


p2v' 


{s-x){q-A)  {s  +  x){q  +  ^) 

(1-c)2 


Using  (37)  and  replacing  corresponding  parts  in  the  second  term  by  the  results  obtained 
above,  we  have  the  following  expression  for  the  derivand  of  the  second  term: 
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2  2 

U  j  m2 


,  ,  pu 

2  2  +  ~  Pi^q i«i  +  a  -  /^2  =  pw  v  +  — 


fa  +  ^)fa  +  j)  fa  “  ^)(s  -  j) 


1  —  c 


+ 


1  -  2c  /  A2  pu 
+ - (yOV'J  +  — 

2  V  4 


(s-jX^-^O  (.s'  +  j)fa  +  A) 


1  -c 


-c(l-c) 


pV 


(.S'  -  x)(q  -  a)  +  fa  +  ^fa  +  'O 


(1-c)2 


1  -2c 


(pv')2  -c(l-c) 


.2.  f 


P  v 


(»-2r)(?-^)  (s  +  dfe  +  2) 
(1-c)2 


+  jUx  -  jU2  =  puv  +  jU{  -  jU2  + 
=  puv  +  Ec  . 


Thus,  a  conservation  law  for  the  chemical  flux  difference  is: 
dp  v  d(puv  +  A) 


dt 


dx 


=  ip  =  Nl  -  N, 


2  ’ 


(51) 


where  lF  is  calculated  as  a  difference  of  the  right-hand  sides  in  (15),  using  the  parametric 
representations  of  the  material-specific  parameters  (16),  (31),  and  (34). 

From  (30)  the  total  thermal  impulse 

p'l  C]  1  +  p' 2  Cj2=  p  Cj 

and  after  summation  of  the  thermal  impulse  equations  in  (9)  and  (10)  we  have  the 
following  expression  for  derivand  of  the  convective  term  with  help  of  (16),  (31),  and  (44): 


p[<hux+  p'2q2u2  +Tx+T2=  p 


fa  +  p v'(l  ~  c))  +  t  —  fa  -  py'c) 


+ 


J2,  .t 


+  2 T  -  ~~  [fa  +  2)(l  -  c)  -  (q  -  X)c\  =  pqu  +  2 T  . 


Thus,  the  conservation  law  for  the  thermal  impulse  takes  the  following  form 

dpq  |  d{puq  +  2T)_  p 

dt  dx 


(52) 


where 


F  =  P- 


^q  +  A  +  q  —  A 

V  ri  r2  7 


=  Pi  +  (Pi 


and  the  arguments  of  the  thermal  relaxation  times  n  and  ri  are  calculated  with  the  use  of 
the  phase's  parameter  representations  (17)  and  (34).  It  is  interesting  to  note  use  of  2  as  a 
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denominator  for  the  flux  of  the  generalised  temperature  in  the  conservation  law  (52)  (for 
an  n-phase  mixture  the  denominator  would  be  n,  i.e.  is  the  number  of  species  in  the 
mixture).  The  occurrence  of  this  denominator  in  the  present  consideration  results  from  an 
association  of  temperature  and  entropy  flux  with  the  total  volume  in  contrast  to  partial 
characteristics  such  as  density  or  pressure. 


The  conservation  law  for  parameter  A  can  be  obtained  by  subtraction  of  the  corresponding 
equation  for  the  thermal  impulse  in  (10)  from  that  in  (9).  The  following  relationship  is  used 


p\qxux  —  T]  —  T2  —  p 


(v  +  ^  (u  +  pV'(±  _  c))_  fa  -  py'c) 


+ 


2.  9 


+  2 Ez  -  ^  ^  [(g  +  /l)(l-c)+(g-/l)c]  =  puA  +  2Ez  . 


Hence,  the  conservation  law  takes  the  form 

dp/ 1  +  d(puA  +  2  A)  _ 
dt  dx 


(53) 


where  A  =  Ex  and 


q  +  A 
T, 


q  -  A 


=  <Pi 


(Pi  ■ 


Similarly,  conservation  laws  for  entropy  s  and  parameter  j  could  be  obtained  from  (14) 
and  (45): 


dpl+d(ml2j) 
dt  dx 


(54) 


and 


dpx  ,  d{pux  +  2B)  _  R  _R  =H 

dt  dx  125 


(55) 


where  B  =  Ex.  The  term  A  was  added  up  in  order  to  take  into  account  dissipation  that  arises 
during  compaction  of  components  of  the  multi-phase  mixture. 

The  energy  conservation  law  could  be  derived  directly  from  the  conservation  laws  (47)- 
(55),  the  equation  of  state  (35)-(36),  and  relationships  for  the  thermodynamic  dependent 
variables  obtained  in  the  previous  Section  (in  fact,  this  procedure  will  be  followed  in  the 
subsequent  Section).  Nevertheless,  to  confirm  the  consistency  of  the  derivations  we  will 
derive  the  energy  conservation  law  for  the  mixture  directly  from  the  phase  conservation 
laws  in  (9)-(10).  To  derive  it,  we  first  obtain  an  auxiliary  relationship: 
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C  (Ml)2  +  (1  -  C)  (M2)2  -  U2  +  p2  C(  1  -  C)(v')2  . 

Using  it,  we  can  evaluate  the  following  part  of  the  energy  flux: 


(  2  x  (  2\ 

t  ^1  f  ^1  f  f 

P\U\  U  +—  +A«  1  U+  —  +PlU\+P2U2  = 

^  ^  i 


=  p  c(u  +  pv'(l-c))  e\+~Y  +(l  -c)(u-pv'c)  e2  +  l-^~  + 


+  {u  +  pv'{\  - c))p\  +  (m  - pv'c)p'2  =  pue  +  p2v'c(l-c)(e1  —e2)  + 

+  —  (cuf  +  (l  - c)u22  )+P  V  — —  (wf  -ul)+up  + p2v'c(l-c)  = 


+  Jy(cwf  +(1-c)m22) 


P\  P2 


=  pu  E  +  —  +/?c(l -c)(/?u')2  m  +  ^ — —  |  +  p2v'c(\ -c)  e{-e2+  —  -  — 

l  2  J  l  2  J  {  Pi  P2 


Using  (44)  and  (45)  along  with  definition  (46),  we  can  expand  the  expression  /1T1  +  J2T2  and 
evaluate  the  total  energy  flux  as  follows 


pu  E  +  —  +Q(//j  -  p2)+ pu  +  pc{\- c)(pv')2  u  +  ^ — —  + 

l  2  )  V  2  ) 

^  Cl2  ( s  +  y  q  +  A  s-r  q  -  A^ \ 

+  2  TJ  +  2AB - - —  •  - + —  •  -  . 


4vc  c  1-cl-c 


Using  the  definitions  of  chemical  potential  (37)  and  generalised  pressure  (38),  we  have  the 
following  compact  expression  for  the  energy  flux 


pu  E  -\ - +  fTA  +  Pu  +  271/  +  2AB  . 

2 


Thus,  the  energy  conservation  law  takes  the  following  form  for  the  two-phase  mixture: 

dp[E  + 11 2  /  )  d  pii\E  +  w 2  /  )+  Pm  +  QA  +  Pu  +  2TJ  +  2  AB 

— - tJd  +  — — 2 - - -J  =  0  .  (56) 

dt  dx 

Summarizing,  the  complete  system  of  equations  for  a  two-phase  medium  involves  the 
mass,  momentum,  and  energy  conservation  laws 


17 


DSTO-TR-1862 


^  +  ^  =  o, 

dt  ox 

dpu  +  d{pu2  +  P)  _ 
dt  ox 


(57) 


dplA  + 


u 


d 


pu  \E  +  uy2  )+  Pu  +  AQ  +  2 TJ  +  2AB 


dt 


dx 


-  =  0 


and  conservation  laws  for  mass  and  volume  concentrations,  the  thermal  impulse,  and  the 
fluctuations  of  the  thermal  impulse,  entropy  and  chemical  flux: 


dpc  d(puc  +  Q) 
dt  dx 

dpO  dpu  6 


0  9 


-  +  - 


=  o 


dt  dx 
dpq  d(puq  +  2  T) 


-  +  ■ 


dt  dx 

dp  A  d(puA  +  2  A)  _ 
dt  dx 

dPA  |  d(puz  +  2B) 
dt  dx 

dpv  d(puv  +  A)  _ 
dt  dx 


=  ~F, 
-G 
=  H , 


(58) 


=  'F  . 


The  system  (57)-(58)  is  closed  by  the  equation  of  state  taken  from  (35)-(36)  in  the  following 
form 


E  =  E(p,c,  6,  q,  A,  s,  j,  v)  .  (59) 

Therefore,  the  thermodynamic  dependent  variables  can  be  calculated  from  the  following 
generalised  thermodynamic  identity 

T  ds-  dE  -  P  dp/p2  -  A  dc  -  A  dj  -  B  dA  -  }  dq  -Q*  dv  -Tldd  .  (60) 

Specifically,  from  (59)  and  (60)  the  closure  relations  for  the  dependent  variables  take  the 
following  form: 

P  -  p2Ep ,  A  -  Ec ,  Tl-Ee,  T  =  Es ,  A-Ex,  J  -  Eq ,  B  =  E\,  El-Ev.  (61) 

The  right-hand  sides  of  the  system  (58)  are  determined  by  the  internal  mechanisms  of  the 
behaviour  of  the  mixture. 

It  is  seen  that  the  jump  conditions  can  be  easily  obtained  for  the  model  because  the 
conservation  laws  (57)-(58)  can  be  used  directly  for  this  purpose.  In  particular,  from  the 
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jump  conditions  it  follows  that  the  dependent  thermodynamic  parameters  P,  A,  T,  A,  J,  B, 
and  Q  and  the  averaged  velocity  u  are  the  parameters  which  are  conserved  through  the 
contact  jump  (when  ignoring  the  parameter  evolution,  associated  with  the  interphase 
friction  and  internal  relaxation  determined  by  the  right-hand  sides  of  (58)). 


6.  Correctness  of  the  model 


The  energy  conservation  law  can  also  be  obtained  from  the  system  (58)-(59)  in  which  the 
energy  conservation  law  has  been  replaced  by  the  conservation  law  for  entropy  (54).  Using 
this  approach  for  the  derivation,  equations  of  the  system  are  multiplied  by  appropriate 
coefficients  and  after  summation  they  result  in  the  energy  conservation  law.  We  will 
conduct  a  backward  equivalent  procedure,  expanding  the  energy  conservation  law  after 
which  the  entropy  conservation  law  is  obtained  and  conditions  of  correctness  can  be 
found. 

Thus,  expanding  the  equation  (56)  and  using  the  mass  conservation  law,  we  can  reduce  it 
to  the  following 

dE  du  dPu  SO  A  dTJ  SAB 

p  —  +  pu  —  + - + - +  2 - +  2 - =  0  , 

dt  dt  dx  dx  dx  8x 

here  d/  dt  =  d/dt  +  u  5/  dx  is  the  particle  derivative.  Expanding  this  formula  with  the  use  of 
the  mass  and  momentum  conservation  laws,  equations  (58),  and  the  relationships  (61),  we 
have 


O 


r  ,  sa^ 

r  so ' 

(  n dT > 

(  „  ^dA\ 

'F - 

+  A 

m(, - 

+  no  +  j 

—  F  —  2  — 

+  B 

-G- 2  — 

k  fa) 

l  dx ) 

l  dx  J 

l  dx  J 

,  rp  dS 

+  pT - b 

dt 


+  A 


H -2  — 
dx 


+  - 


SOA 

dx 


„  dTJ  „  SAB 
+  2 - +  2 - =  0. 


dx 


dx 


Replacing  the  particle  derivative  of  entropy  from  (54)  and  observing  cross-annihilation  of 
space  derivatives,  we  can  eventually  obtain  a  restriction  on  the  right-hand  sides  of  the 
system  (58)  and  equation  (54): 

QT  +  Am0  +  no  -JF-  BG  +  T(Rx  +R2+A)+AH  =  0  . 


Now,  we  can  rewrite  denotations  of  Ni  and  Ni  in  (15),  when  using  T1  =  N\  -  N2,  through  Ri 
and  R2  in  (14)  as  follows 


N,= 


n0  -uxm0 
pc 


qlRl+qlslm0+<plsl  ,  N2  = 


n0-u2mQ 

T.  \  ^2-^2  ~  q 2S2m0  (PlS2  • 

p{l-c) 
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Expanding  F,  G,  and  H  with  the  substitution  of  their  expressions  in  (52),  (53),  and  (55),  we 
have 

+  A/Wq  +  no  —  +  (p~, )  —  —  (p~, )  +  +  R~,  +  a)  +  A(i?j  —  R~, )  =  0 . 


Using  the  relationships  (37)  and  (46)  for  A  and  Q  and  presenting  the  averaged  variables 
back  through  the  variables  for  each  phase  (material-specific  variables)  with  the  help  of 
equations  (44)  and  (45),  from  the  last  two  formulas  it  follows 


Q 


n()-u,m()  n(.  -  ujn,. 

~  VA  +  V«0  +  <PlS  1  + - U  "  \  +  <hR2  +  <hS 2m0  ~  V2S2 


pc  p{\  -  c) 

(l-2  c\pv')2  (  ^ 

Mi  ~Mi  +~ - J - ^W^i+q2s2) 


+ 


m0  +  no  - 


-  Jl(pl  -  Qsl(pl  -  J2(p2  +  Gls2(p2  +  Jj  A,  +  QRlql  +  T2R2  -  QR2q2  +  JA  =  0  . 


After  obvious  simplifications  we  have 


pv'c{\  -  c) 


n0  -  utm0  n0  -  u2m0 


1  -c 


+  m0(  1  -  2c)(/u/')  +  m0pnv'  -  u]n0  +  u2n0  + 


+  no  +  TA  =  0 


that  is  finally  reduced  to 

no  +  TA  =  0  . 


Thus,  fulfilment  of  the  energy  conservation  law  imposes  the  following  restriction  on  the 
compaction  dissipation  term  introduced  in  (55): 


no  _  pU2 

~T~1\ 


(62) 


Therefore,  as  soon  as  the  exchange  terms  lo,  mo,  and  no  in  (9)-(10)  are  chosen  in  such  a  way 
that  Ri  +  R2  >  0  in  (14),  then  the  total  entropy  production  is  non-negative  from  (62). 

For  an  analysis  of  when  the  system  (57)-(58)  is  hyperbolic,  we  will  conduct  it  in  the  way 
performed  in  [15]  for  the  single-temperature  model.  The  eigenvalue  analysis  is  fairly 
complex  for  the  system  (57)-(58)  that  can  be  written  in  the  following  matrix  form 


dU  AdU 

- +  A - 

dt  Gx 


=  B, 


(63) 
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where  U  T  =  {  p,  c,  9,  u,  v,  q,  A,  s,  j  }.  The  point  of  the  approach  used  in  [15]  was  to 
reformulate  a  final  system  written  in  the  mixture  variables  into  an  equivalent  form, 
employing  material-specific  variables  for  the  purpose  of  the  eigenvalue  analysis.  For  the 
present  case  the  latter  set  of  variables  can  be  assembled  as  the  following  vector: 

U  T  =  {0,  pi,  ui,  qi,  si,  pi,  ui,  qi,  S2}  •  (64) 

The  compaction  dynamics  equation  (49)  can  be  reduced  to  the  following  one: 


86  86  ® 

- h  U - =  —  . 

8t  8x  p 


(65) 


We  remind  that  the  final  system  (57)-(58)  (which  could  be  used  in  an  alternative  form  with 
the  energy  conservation  law  from  (57)  replaced  by  the  entropy  equation  (54))  is  a  result  of 
equivalent  proper  combinations  of  the  systems  (9),  (10),  and  equations  (14)  and  (15).  In 
order  to  manipulate  with  the  'real'  densities  for  each  phase  we  use  the  following 
relationships: 


p\  =  01  pi  ,  p'i=  dip 2  , 


where  denoted  6i  =  6  and  62  =  1  -  6  .  Keeping  in  mind  the  following  relationship  for  the 
averaged  density  p  =  dipi  +  9-ipi,  we  can  derive  from  the  continuity  equation  (47)  and  from 
the  equation  (48)  for  pc  =  pT  =  0ip\  the  following  equations  for  the  phase  densities: 


8p\  |  dpxux  |  a 


(ip  -u)  86 


8t 

dpi 

8t 


8x 

dp2u2 

8x 


6X  8x 


A® 


P  J 


=  M  1 


p2(u2-u)Sd  1 


6 , 


ox  6-, 


-  m0  + 


Pi® 

P  J 


(66) 


=  AT, 


From  the  entropy  conservation  law  (54)  and  the  conservation  law  for  the  entropy 
disequilibrium  j  we  can  obtain  the  material-specific  entropy  conservation  laws  similar  to 
(14)  with  an  additive  of  A/ 2  to  the  right  hand  sides.  From  these  equations  we  can  derive 
auxiliary  equations  for  0iSi  and  62s2  with  the  use  of  relationships  p Tsi  =  pidiSi  and  p'2s2  = 
pidisi  and  equations  (66): 


Pi 


Pi 


86,s<  86,s< 

+  U, 


8t 


8x 


A  /  \86  8JX  A 

-Pi*  Mi  = 

6x  6x  2 


86, s,  86tS, 

— —  +  u2 — — 


8t 


8x 


A  /  \86  8J2  A 

+  P 2S 2 2  U)  ~  n  ~  ^2  +  ~  &2S 2^^  * 

OX  OX  2 


(67) 


Similarly,  from  the  conservation  laws  (52)  and  (53)  for  the  thermal  impulse  q  and  its 
disequilibrium  A,  we  can  obtain  the  conservation  laws  for  the  phase's  thermal  impulses 
that  are  identical  to  the  thermal  impulse  conservation  laws  in  (9)  and  (10).  After  some 
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manipulation  for  exclusion  of  the  volume  concentration  they  give  us  another  pair  of 
auxiliary  equations  for  p\Cf\  and  p2C]2- 


A 


dpxgx  |  dplqiul 
dt 


dx 


^  (  280  dTx 

+  pxq\ux  -u)—  +  —  =  -(px  - 

ox  ox  p 


J 


8p2,;+8A^^_^(^_u)8g+8r2 


dt 


dx 


J 


dx  dx 


~(pi  + 


P 


(68) 


After  multiplication  of  the  first  equation  of  (67)  by  i]i  and  the  first  equation  of  (68)  by  Si 
followed  by  their  summation,  we  can  obtain  an  equation  for  pdhqiS]  and,  similarly,  for 
pidiqiSi. 


SpAqi^  ,  dpxuxexqxsx 


dt 


=  q  i 


dx 


dT :  a/, 

+  5, - 1-  q, - — 

dx  dx 


A 


A 


Rx  +  —  -OxsxMx 

2  j 


+  s. 


Vx 


P  J 


=  H 


1  5 


dp262q2s2  dp2u202q2s2  dT- 


dt 


=  q2 


dx 


+  S7 


dx 


+  q  2 


dJ , 


(69) 


dx 


f 


A 


A 


R2  +  ^  02s2M2 


+  s- 


■cp2  + 


AA° 
P  J 


=  h2. 


Rewriting  the  conservation  law  (51)  for  the  chemical  flux  pv  in  its  equivalent  form  for  the 
material-specific  parameters  (as  a  difference  of  equations  in  (1)),  and  then  adding  it  to  the 
first  equation  from  (69)  followed  by  subtraction  of  the  second  equation  of  (69),  we  can 
obtain  the  equation  for  the  velocity  difference  in  the  following  form: 


d(u j  -  u2 ) 


(  2  / 

UX  /  U 2 


dt 


+  - 


2  A 

2J 


d(p,  -  Pi )  ,  „  dTx 


dx 


dx 


+  5, 


dx 


+  q  i 


dJx 

dx 


-s. 


dT 2_ 
dx 


q2 


dJ2 

dx 


=  XV  +  Hx  -  H2  . 


Expanding  the  derivative  of  the  chemical  potentials  according  to  (6),  this  equation  can  be 
reduced  to  the  following: 


du,  du,  du ,  du1  1  dp, 

— 2 - -  +  u,  — --u,  — -  + - — 

dt  dt  dx  dx  px  dx 


1  dp2 

p2  dx 


=  X¥  +  HX 


-H 


2  • 


(70) 


The  momentum  conservation  law  (50)  can  be  rewritten  in  the  form  of  a  sum  of  the 
momentum  conservation  laws  from  (9)-(10).  Excluding  the  variables  p\  and  dx  from  the 
derivands  of  (50),  we  obtain  the  following  equation: 
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®\P\ 


du  j 
dt 


■  +  M, 


1 


■  + 


dx  dx 


+  02p2 


dip 

dt 


■  +  Ur, 


du , 


■  + 


1  dp2 


dx  p2  dx 


+  {Px = 

dx 


Combining  appropriately  this  equation  with  (70),  we  can  obtain  equations  for  the 
velocities  of  each  phase: 


du,  du,  1  dp, 

— 1  +  u,  — 1  + - —  + 

dt  dx  px  dx 

du2  du2  1  dp2 

b  U  2  I 

dt  ~  dx  p2  dx 


(Px  -Pi) 


do_ 

dx 


+ 


(A  -Pi) 


oo_ 

dx 


02p2^  +  Hx-H2)=Gx, 

=  ~9\P\ +  Hl-H2)=G2 . 


Finally,  from  (67)  and  (68)  we  can  derive  equations  for  q,  and  s,: 


dqx  dqx  1  8TX 


dt 

Sq2 

dt 

dsl 

dt 

dp 

dt 


+  u 


1 


dx  px0x  dx  px9x 


<Px 


PxO 


qxM  i 


+  u 


dq2 


+  ■ 


1  8T ; 


1 


dx  p262  dx  p202 


+  u. 


'  +  li  0 


dsx 

dx 

dSr, 


+  ■ 


1  dJx 
px6x  dx 
1  8J 


1 


■<Pi  + 


A 


P  ) 

PiOi® 

P  ) 


Px 


=  Fx 


q2M2 

Pi 


=  F- 


2  ’ 


+  - 


pA  v 

l 


Rx  +  --6lslM  j 
2  7 


1  > 


7 


dx  p202  dx  p20 


2  V 


i?2  +  2  @2SlM  2 


+  - 


s2® 

pA 


=  AC 


(71) 


(72) 


Summarising,  from  (65),  (66),  (71),  and  (72)  we  can  obtain  a  system  of  equations,  which  is 
equivalent  to  the  system  (57)-(58),  in  the  following  form: 


d9  dO  ® 

—  +  u  —  =  —  , 
dt  dx  p 


dpx  dpx  dux  dO 


dt 

dux 

dt 


■  +  u. 


■  +  u. 


dx 

du. 


+  Px^  +  a  x~^~Mx 
Ox  dx 


1  dp 
dx  px  dx 


■  + 


'+6^  =  0, 
dx  1 


A.+„iAl+_LAI  =  F| 

dt  dx  P\0X  dx 


dsx  Os,  1  dJ, 

— -  +  u,  — L  + - ~  =  NX 

dt  dx  px0x  dx 


(73) 


23 


DSTO-TR-1862 


■  +  U-, 


dpi  ,  ..  dp2 
dt 
du  i-i 


dt 


-  +  un 


dx 

du. 


du 2  dO 

+  P2—^-a2  —  =  M2  , 

dx  dx 


1  dp  2 

dx  p2  dx 


■  + 


+  b^  =  G 


dx 


2  ■> 


dq2  dq2 


dt 


■  +  M, 


+  - 


1  dT ; 


ds2 

3? 


■  +  W'j 


dx  p202  dx 

ds,  1  dJ, 


=  F. 


2  ’ 


+  - 


dx  p262  dx 


=  N. 


2  ’ 


where  ai  =  pi(ui  -  u)/6i  and  ca  =  pi{u2  -  u)/ 02,  and  b  =  pi  -  pi.  We  assume  that  the  equation 
of  state  for  each  phase  has  independent  mechanical  and  thermal  impulse  parts  as 
proposed  in  [12]: 

e\  (A  -  e\  (a  j  (<?i  )  j  A  (A  d2>?i)_  ^2  (  A  ’A  )  “*~  A  (#2  )  • 

We  will  use  the  following  denotations  (specified  just  for  the  first  phase): 
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(74) 


Thus,  a  presentation  of  the  equivalent  system  (73)  in  the  matrix  form  (63)  for  the  vector  of 
variables  (64)  can  be  obtained,  where  the  matrix  A  is: 
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and  the  right  hand  side  in  (63)  is 
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B  1  -  {O/ p,  Mi,  Gi,  F i,  Ni,  M2,  Gi,  F2,  N2} . 

Eigenvalues  t]  of  this  system  can  be  found  from  the  following  characteristic  equation 

det(A  -  77/)  =  0  .  (76) 

Due  to  the  structure  of  the  matrix  A  (75),  the  characteristic  equation  (76)  can  be 
decomposed  into  the  following  product: 
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fo2-cr). 


£1- 


Q.-s,  rfe 

{pAY J  &Y 


fe-Q2)- 
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02»2  ^ 


(pA)2)  (02) 


1 2  ^2 


=  0,  (77) 


where  £0  =  u  -  fj,  £1  =  U\  -  rj,  and  £2  =  u2  -  rj.  The  solution  £0  =  0  corresponds  to  the  evolution 
of  the  volume  concentration  and  two  multipliers  in  (77)  correspond  to  the  characteristic 
equations  for  the  first  and  second  phases.  Analysing  one  of  them,  the  conditions  of  the 
existence  of  the  roots  are 


(S)  .  77  f2  . 77 
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Recalling  denotations  (74),  we  can  rewrite  the  conditions  as  follows 
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After  introduction  of  the  specific  volume  as  V\  =  1/ pi,  the  conditions  take  the  following 
well-known  form 


-j2  T  ^2  X  rs2„X 

d  ei  ^  Q  d  ex  d  ex  > 


dq[ 


8V,  ds 


(78) 


which  represents  the  traditional  conditions  of  the  equation  of  state  convexity  for  e\ x  and 
eiT.  The  same  restriction  should  obviously  be  imposed  for  e2x  and  e2T  as  well. 


7.  Discussion  and  conclusions 


The  system  of  equations  (57)-(58)  has  been  obtained  in  a  compact  form  for  a  two-phase 
heat-conducting  medium.  It  appears  that  a  generalisation  for  the  multi-phase  case  can  be 
conducted  in  an  obvious  way  similar  to  the  generalisation  conducted  in  [2] .  The  conditions 
of  hyperbolicity  of  the  system  (78)  have  been  reduced  to  the  requirement  of  convexity  for 
the  mechanical  and  thermal  parts  of  the  equation  of  state  for  individual  phases.  It  should 
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be  noted  that  non-strict  hyperbolicity  might  be  associated  with  the  phase  degeneration  and 
a  degeneration  of  the  heat  conductivity  (which  may  happen  when  Hi=0  in  (74))  that  is  out 
of  the  scope  of  the  report.  The  correctness  of  the  model  has  been  linked  with  the 
requirement  of  non-negative  dissipation  (62)  associated  with  the  compaction  dynamics.  It 
is  interesting  to  note  that  the  structure  of  the  characteristic  equation  (77)  preserves  the 
characteristic  structure  of  the  original  single-phase  model  when  expanding  it  to  the  two- 
phase  case.  This  allows  one  to  solve  a  Riemann  problem  for  a  two-phase  mixture  with  the 
help  of  single-phase-model  Riemann  solvers  for  each  phase,  facilitating  development  of  a 
suitable  numerical  algorithm  based  on  the  Godunov  method. 

It  should  be  noted  that  the  generalisation  of  the  model  to  the  multi-dimensional  case  does 
not  look  obvious,  and  the  major  reason  is  the  difficulties  in  a  tensorial  generalisation  of  the 
equation  of  state. 

The  appearance  of  the  denominator  2  (that  presents  the  number  of  species  in  a  two-phase 
mixture),  entering  the  governing  equations  (54)  for  entropy,  equations  in  the  system  (58) 
for  the  thermal  impulse  and  for  fluctuations  of  the  thermal  impulse  and  entropy,  and  the 
energy  conservation  law  in  (57)  can  be  interpreted  by  the  fact  that  the  dependent 
thermodynamic  parameters  T  and  /  are  not  really  partial  parameters,  in  contrast  to  density 
and  pressure. 
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